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Abstract 

Highly optimized tolerance is a model of optimization in engineered systems, which gives rise 
to power-law distributions of failure events in such systems. The archetypal example is the highly 
optimized forest fire model. Here we give an analytic solution for this model which explains the 
origin of the power laws. We also generalize the model to incorporate risk aversion, which results 
in truncation of the tails of the power law so that the probability of disastrously large events is 
dramatically lowered, giving the system more robustness. 



1 



In a series of recent papers, Carlson and Doyle 0, [| have proposed a model for 
designed systems which they call "highly optimized tolerance" or HOT. The fundamental 
idea behind HOT is that systems designed for high performance naturally organize into 
highly structured, statistically unlikely states that are robust to perturbations they were 
designed to handle, yet fragile to rare perturbations and design flaws. As an example they 
consider an idealized model of forest fires [[J . In this model a forester is charged with finding 
the optimal distribution of the trees on a grid so as to maximize tree harvest in the face of 
occasional fires that burn complete connected clusters of trees and are started by sparks that 
arrive with a given spatial distribution. They find that optimizing the harvest, or yield, for 
the model gives rise to a segmented forest consisting of contiguous patches of trees separated 
by firebreaks, and that the resulting distribution of fire sizes usually follows a power law. 
While this type of configuration does typically achieve very good yields, the system is also 
fragile in the sense that perturbations to the firebreaks or changes in the spark distribution 
can lead to substantially sub-optimal performance. They argue that these are pervasive 
phenomena: high-performance engineering leads to systems that are robust to stresses for 
which they were designed but fragile to errors or unforeseen events. 

In this paper we argue that simple yield maximization is problematic even if there are no 
errors in firebreaks or changes in the spark distribution. Because the power-law distributions 
generated by yield maximization have fat tails, disastrously large forest fires occur with non- 
negligible frequency — far greater frequency than one would expect from intuition based on 
normal distributions. This idea, that yield optimization can lead to ruinous outcomes, is not 
new. For the classic problem of gambler's ruin, for example, it is well known that optimizing 
total return leads to ruin with probability one. By contrast, if one is willing to accept 
suboptimal returns it is possible to construct gambling strategies that are immune to ruin [|J . 
Applying similar ideas in the present context, we show that a risk- averse engineer who is 
willing to accept some loss in average system performance can effectively limit the large 
deviations in the event size distribution so that disasters are rare. We call this variation on 
the HOT theme "constrained optimization with limited deviations", or COLD. By avoiding 
total ruin, a COLD design is more robust than a HOT one, even in a world of perfect 
error-free optimization. 

To demonstrate the difference between HOT and COLD we first revisit the HOT forest 
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fire model. We give an analytic solution for the model which shows that the distribution of 
fire sizes does indeed follow a power law, the cumulative distribution having a exponent — (1 + 
1/d), where d is the dimensionality of the system. Using a fast percolation algorithm || we 
perform numerical simulations that confirm the value of this exponent. We then generalize 
our solution to include risk aversion in the design, thereby breaking the power law scaling 
and dramatically reducing the frequency of disastrously large events. 

Following Refs. [I] and § then, we consider a forest divided into a large number of regions 
or patches, with firebreaks between them that prevent the spread of fire from one patch 
to another. Although the original forest fire model was based on a lattice, the model we 
consider is a continuum one, since this makes the mathematical treatment more tractable. 
For large system sizes, we expect the behavior of this continuum model to converge to that 
of the lattice model. 

It is assumed that during the lifetime of the forest a single spark lands at a random 
position r and starts a fire that burns the surrounding patch. Let us denote the area of this 
patch by s(r). The forest is then harvested, giving a yield equal to the area of the remaining 
forest. In units where the total area of the forest is one, the yield is 1 — s(r) — F, where F 
is the cost in terms of yield of constructing the firebreaks. 

Because dimensionality is an important property of HOT systems, we consider the model 
for general dimension d. If the cost of constructing firebreaks is a per unit length (or per unit 
surface area for d > 2), then the cost of the firebreak surrounding a patch m is agdsm~ 1 ^ d , 
where s m is the value of s(r) in patch m and g is a geometric factor of order 1 that depends on 
the geometry of the lattice and the shape of the patches. In the lattice version of the forest 
fire model, a is simply equal to the lattice parameter (i.e., the nearest-neighbor spacing), 
but in the continuum model we are at liberty to give a any value we feel to be appropriate. 
As we will see, as long as a is finite its value does not affect the shape of the distribution of 
fire sizes. 

Because s(r) is constant inside each patch, the integral of l/s(r) over any patch is iden- 
tically 1, and hence, summing over all patches, the total area occupied by firebreaks is 

F = agdJ2s%- 1)/d = agdJ2^~ 1)/d [ TT = W dV (1) 
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Letting the normalized probability distribution of sparks be p(r) , the mean yield is then 

Y = 1 - J p(r)s(r) d d r - agd J s(r) -1/d d\ (2) 

where the integrals run over the entire area of the forest. 

To find the maximum yield with respect to the patch sizes s(r), we set the functional 
derivative 5Y/5s(r) = 0, giving 

r -| d/(d+l) 

f \ a 9 /oN 

_p(r)_ 

for all r. The optimal yield is then given by substituting back into Eq. (0) to get 

Y opt = l-(d + l){ag) d '^ j p(r) 1 /^ d d r, (4) 



and the optimal number of patches is 



(ag)- d ^ / p(r)^ d+1 ) d d r. (5) 
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For the lattice version of the model, a goes as L^ 1 , where L is the (linear) system size, and 
hence the number of patches should scale as L d ^ d+1 \ i.e., as L 1 / 2 in one dimension or L 2//3 in 
two. Numerical experiments on a one-dimensional system confirm this, giving a measured 
exponent of 0.47 ± 0.03. 

Now we wish to calculate the distribution p(s) of fire sizes that arises if we make this 
choice of patch sizes. We have 

/ \ , A d r d d rdp d + 1 d d r _( 2 +i/<fi ^ 

p(s) = p(r) — = p(r) — - = -ag—ptf—s (6) 

where d d r here represents the volume of the space between the contours p and p + dp on the 
p(r) surface. As we will show, the term p(r) d d r/dp is constant or contributes logarithmic 
corrections for a wide selection of possible distributions p(r) , while the principal power-law 
behavior in the event size distribution comes from the factor s~( 2+1 / rf ) ||. 

An alternative method for deriving Eq. (|^) is to maximize the simple yield functional 
Y = 1 — j p(r)s(r) d d r, subject to a constraint that fixes the volume F occupied by the 
firebreaks (Eq. (|T])). This method, which is similar to the approach taken in Ref. [I], is 
equivalent to the method above, via a Lagrange transform, provided F is chosen so as to 
make the corresponding Lagrange multiplier equal to agd. 



Consider then the case (which covers all the examples in Refs. [I] and |2|) of a distribution 
of sparks with a single maximum at the origin, so that the volume d d r takes the form of an 
annulus enclosing the origin. If we denote by Qd a (i-dimensional solid angle centered on the 
origin, then a volume element of the annulus is r dr dQd- In terms of the thickness dp of 
the annulus, we can write dr = (r dp) / (r ■ Vp), and integrating our volume element over the 
contour p = constant we find 

d d r f r d d£l 
dp 



(7) 
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For example Carlson and Doyle studied the case of a spark distribution in two dimen- 
sions having the form of the product of two Gaussians with different widths: 



p(r) = iVexp 
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where N is a normalization constant. For this distribution the denominator of the integrand 
in Eq. fl?]) is 
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which is constant over our contour of constant p. The element of solid angle in two di- 
mensions is simply the element of polar angle d8, and hence Eq. ([?]) simplifies in this case 
to 



d 2 r 
dp 
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where A(p) is the area enclosed by the contour. This contour is a line of constant (xl<j x ) 2 + 
(y/a y Y, i-e., an ellipse, which has major and minor axes a = \/2a 2 \og(N/p) and b = 
\f2a y 2 \og(N/p). Thus the area enclosed by the contour is A{p) = nab = 2na x a y \og{N/p). 
Combining Eqs. (|Bp and (0) we then find that the distribution of event sizes is 



p(s) = 37Tcr x (Xyag s 5 ^ 2 . 



(11) 



Thus, for the Gaussian case in two dimensions the model generates a perfect power-law with 
slope — |. 
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FIG. 1: Cumulative distribution of fire sizes for simulations of the forest fire model in one and two 
dimensions, plotted on logarithmic scales, (a) One dimension with size 10000 and an exponential 
spark distribution, (b) Two dimensions with size 128 x 128 and a Gaussian spark distribution. 
The results have been averaged over a number of different values of the parameters of the spark 
distributions to improve the statistics. The dotted lines show the expected slopes of —2 and — |. 

This argument is easily generalized to other spark distributions and other dimensions. 
We find that the HOT forest fire model generates a perfect power law with slope —(2 + 1/d) 
for all dimensions when a spark distribution of the form p(r) = iVexp(— Y2t=i[ x i/' 7 i\ d ) * s 
used. As a test of this prediction, we show in Fig. [I] numerical results from direct simulations 
of the forest fire model in one and two dimensions with distributions of this type |§. For 
better visualization and analysis the distributions pictured are cumulative, so the expected 
slope is —(1 + l/d), rather than —(2 + 1/d). As the figure shows, the slopes of the observed 
distributions are in good agreement with this prediction. 

We note in passing that the slope of — (1 + l/d) for the cumulative distribution of fire 
sizes seen in both our exact solution and our numerical results is different from the slope 
of approximately —1 found numerically by Carlson and Doyle in two dimensions. The 



source of this discrepancy is unclear, although it may be that the simulations of Ref. |2| 
provided too few data points to make an accurate evaluation of the exponent possible. We 
note also that the value — | for the two-dimensional case is quite different from the slope 
of — | measured for the cumulative size distribution of real forest fires ||, |j. 

Other functions with exponential tails also generate power laws, but give logarithmic 
corrections as well. For instance, if p(r) = iVexp(— Ylt=il x i/ a iV) w hh 7 ^ d then Eq. ([T0|) 
still applies, but now A(j>) ~ [\og(N / p)] d ^ and hence 

p(r)^ ~ [\og{N/ P )} d ^-\ (12) 

Thus the distribution of event sizes fundamentally still follows a power law with slope — (2 + 
but there is a logarithmic correction. Similar logarithmic corrections are noted in 
Ref. 0. 

Spark distributions with power-law tails also (unsurprisingly) give power-law event size 
distributions, but in this case the exponent of the distribution is non-universal, varying with 
the exponent of the spark distribution. For example, if p(r) takes the generalized Lorentzian 
form 

N 

pir) = E,teM)" + r- (13) 



1 YA , (14) 



then we find 



dp p(r) \ p 

which goes asymptotically as p~ d l v in the power-law tail where p(r) becomes small. Thus 
the tail of the distribution of event sizes goes as p(s) ~ s -(2+i/d-d/i/) _ 

We now turn to the COLD variant of the forest fire model, which incorporates risk 
aversion. In constructing this model we are guided by theories of risk aversion in economics, 
where the subjective benefit of outcomes is typically a nonlinear function of the loss s, 
which is captured by a utility function u(s) 0. Sensible utility functions are decreasing 
with increasing loss: u' < ||10|| . Risk aversion also implies that u" < 0, so that the negative 
utility of bad outcomes is weighted more strongly than the positive utility of good outcomes. 
One standard family of utility curves that achieves this is the one-parameter family 

u(s) = (15) 
a 
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Note that, since we will be concerned only with maximizing utility, u(s) is arbitrary to 
within both additive and multiplicative constants. For a = 1, Eq. flT5| ) gives u = 1 — s and 
maximizing utility is precisely equivalent to minimizing loss. For a < 1, we have risk-averse 
utility functions and for a < we are infinitely averse to losing our entire investment. 
Our goal now is to maximize the average utility functional 

U = J p(r)u{s{r)) d\ (16) 

subject to the constraint of fixed F, Eq. ([!]) (or equivalently maximize a combined utility 
functional similar to Eq. (0)). Carrying out the functional derivatives and using the utility 
function of Eq. (|lT>D, we find that the optimum U corresponds to 

p(r)s(rY d+1)/d [l - sir)}*- 1 = A, (17) 

where A is a Lagrange multiplier whose value can be calculated from Eq. ([I]). The distribution 
of event sizes is given by Eq. (^j) as before, and using Eq. ( |T7D we find that the derivative 
dp/ds, which gives the principal variation in p(s), is 

dp _ (q + l/d)s - (1 + 1/d) 

ds (1 - s) a s 2 +V d ' 1 J 

For a = 1 our utility maximization is equivalent to simple yield maximization (HOT), so 
it is not surprising to observe that when we set a = 1 in the above expression we recover 
our previous s~( 2+1 / d ) power-law. For a < 1, we have risk-averse utility functions (COLD), 
which give rise to event distributions following the s~( 2+1 / d ) form for small event sizes, but 
having lower probability of large event sizes. When a < 0, event probability tends to zero 
as s — > 1, as we would expect. 

In Fig. ^ we compare the distribution of event sizes in HOT and COLD regimes for 
a variety of values of the risk- aversion parameter a. The figure shows that the COLD 
distribution approaches the HOT one as a approaches 1. For a large and negative the HOT 
power law is followed for only a small portion of the range of event sizes — about 20% in the 
case of a = —5. 

It is worth noting that while risk aversion truncates the power-law behavior in the event 
size distribution, the distribution of the utilities of events still follows a power law: we find 



S 



0.0 0.2 0.4 0.6 0.8 1.0 

fractional size of event 



FIG. 2: Event size distributions for HOT and COLD regimes in two dimensions. The dotted line 
is the distribution for the HOT regime (a = 1) and the solid lines are for the COLD regime with 
(top to bottom) a = —1 ... — 5. The distributions are not normalized — they cannot be since they 
diverge at the origin. In practice Eq. (jl]) provides a lower cutoff on s and makes the distribution 
normalizable. Inset: the same data on log-log scales. 

that the tail of the distribution of utilities goes as p{u) ~ u^ 13 with /3 = (2a — l)/a. Note 



that this exponent is independent of the system dimension d [|Tl[ . 

While the introduction of the utility function has reduced the risk of large losses, the 
optimal utility solution does not normally coincide with the optimal yield solution, and 
hence we pay a cost for risk aversion in terms of yield. For the lattice forest fire model 



however we find that the cost paid is small [12|. For example, in a 128 x 128 two-dimensional 



system with a Gaussian spark distribution, we find numerically that the mean yield at the 
a = 1 optimum (HOT) is 0.904, dropping to 0.900 for a = -3 and 0.888 for a = -5. It 
appears therefore that the introduction of risk aversion garners substantial benefits in terms 
of the reduction of large losses — and the complete elimination of 100% losses — while at the 
same time costing us only a few percent at most in terms of average system yield. 
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We conjecture that the suppression of power law tails in the COLD event size distribution 
will also make the system more robust against the other problems mentioned in the intro- 
duction, namely errors in the design and changes in the spark distribution. The truncation 
of the power law means that the largest patches in the COLD solution are considerably 
smaller than those in the HOT solution. Thus, if a design flaw, such as a gap in one of 
the firebreaks, causes two patches to merge, the resulting combined patch is smaller too. 
Similarly, if the spark distribution is changed, the size of the resulting fires is smaller, and 
hence the effect on the average yield is not as catastrophic as in the HOT case. (A similar 
conjecture is also made in Ref. [j].) 

In this paper we have examined in detail only the forest fire model, but similar principles 
should apply to other problems as well. We have shown that in order to produce power-law 
event size distribution, the HOT model requires the auxilliary assumption of risk-neutrality. 
If humans are risk-averse they will tend to prefer COLD designs, although this does not 
necessarily mean that HOT designs never occur. It might be, for instance, that blind 
evolutionary processes of the type found in natural systems would simply optimize yield, 
without risk aversion. On the other hand, COLD designs are more robust to rare events 
than HOT designs, and therefore might be selected for on long time-scales. Of course, in 
the real world, imperfect designs that fail to optimize either yield or utility are always a 
possibility too. 
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